module response

use, intrinsic :: iso_c_binding
use basic
use avd
implicit none
integer, parameter :: MPR = 30

contains

subroutine newmark(acc,n,dt0,zeta,P,SPA,SPI,SPV,SPD,SPE) bind(c)
!     Newmark-β 法求解 SDOF 的峰值响应
    integer, intent(in) :: n
    real(8), intent(in) :: dt0, zeta
    real(8), intent(in) :: acc(n), P
    real(8), intent(out) :: SPA, SPV, SPD, SPE
    integer, intent(out) :: SPI

    real(8) :: rc(3),rl(3),ac,al,aa,ein,da,rvl
    real(8) :: beta, gamma
    real(8) :: k,c,w,keff,kinv,feff,dt
    real(8) :: b1, b2, b3, b4, b5, b6, b7, b8
    integer :: i,j,r

    beta = 1.d0/4.d0
    gamma = 0.5d0

    if ( dt0*MPR > P ) then
        r  = ceiling(MPR*dt0/P)
        dt = dt0/r
    else
        r = 1
        dt = dt0
    end if

    b1 = 1.d0/(beta*dt*dt)
    b2 = 1.d0/(beta*dt)
    b3 = 1.d0/(beta*2.0d0)-1.d0
    b4 = gamma/(beta*dt)
    b5 = gamma/beta-1.d0
    b6 = 0.5d0*dt*(gamma/beta-2.0d0)
    b7 = dt*(1.0d0-gamma)
    b8 = dt*gamma

    w = pi2/P
    k = w*w
    c = 2.0d0*zeta*w
    keff = k+b1+b4*c
    kinv = 1.0d0/keff

    !  initiation
    rl = 0.d0
    rc = 0.d0
    SPA = 0.0d0
    SPV = 0.0d0
    SPD = 0.0d0
    SPI = 0
    ein = 0.0d0
    al = 0.d0
    rvl = 0.d0

    !  computation of response
    do i = 1, n, 1

        ac = acc(i)
        da = (ac - al)/dble(r)

        do j = 1, r, 1
            ac = al + da*j
            feff = ac+(b1*rl(1)+b2*rl(2)+b3*rl(3))&
            +c*(b4*rl(1)+b5*rl(2)+b6*rl(3))
            rc(1) = feff*kinv
            rc(2) = b4*(rc(1)-rl(1))-b5*rl(2)-b6*rl(3)
            rc(3) = ac-k*rc(1)-c*rc(2)
            rl(1) = rc(1)
            rl(2) = rc(2)
            rl(3) = rc(3)
        end do

        ein = ein+0.5d0*dt0*(rvl*al+rc(2)*ac)
        aa = -rc(3)+ac

        if (abs(aa) > abs(SPA)) then
            SPA=aa
            SPI=i
        endif

        if (abs(rc(2)) > abs(SPV)) then
            SPV=rc(2)
        endif

        if (abs(rc(1)) > abs(SPD)) then
            SPD=rc(1)
        endif

        al = ac
        rvl = rc(2)

    end do
    SPE = sqrt(2.d0*ein)
    return

end subroutine newmark

subroutine rnmk(acc,n,dt0,zeta,P,ra,rv,rd) bind(c)

    integer, intent(in) :: n
    real(8), intent(in) :: dt0, zeta
    real(8), intent(in) :: acc(n), P
    real(8), intent(out) :: ra(n),rv(n),rd(n)

    real(8) :: rc(3),rl(3),ac,al,da
    real(8) :: beta, gamma
    real(8) :: k,c,w,keff,kinv,feff,dt
    real(8) :: b1, b2, b3, b4, b5, b6, b7, b8
    integer :: i,j,r

    beta = 1.d0/4.d0
    gamma = 0.5d0

    if ( dt0*MPR > P ) then
        r  = ceiling(MPR*dt0/P)
        dt = dt0/r
    else
        r = 1
        dt = dt0
    end if

    b1 = 1.d0/(beta*dt*dt)
    b2 = 1.d0/(beta*dt)
    b3 = 1.d0/(beta*2.0d0)-1.d0
    b4 = gamma/(beta*dt)
    b5 = gamma/beta-1.d0
    b6 = 0.5d0*dt*(gamma/beta-2.0d0)
    b7 = dt*(1.0d0-gamma)
    b8 = dt*gamma

    w = pi2/P
    k = w*w
    c = 2.0d0*zeta*w
    keff = k+b1+b4*c
    kinv = 1.0d0/keff

    !  initiation
    rl = 0.d0
    rc = 0.d0
    ra = 0.d0
    rv = 0.d0
    rd = 0.d0
    al = 0.d0

    !  computation of response
    do i = 1, n, 1

        ac = acc(i)
        da = (ac - al)/dble(r)

        do j = 1, r, 1
            ac = al + da*j
            feff = ac+(b1*rl(1)+b2*rl(2)+b3*rl(3))&
            +c*(b4*rl(1)+b5*rl(2)+b6*rl(3))
            rc(1) = feff*kinv
            rc(2) = b4*(rc(1)-rl(1))-b5*rl(2)-b6*rl(3)
            rc(3) = ac-k*rc(1)-c*rc(2)
            rl(1) = rc(1)
            rl(2) = rc(2)
            rl(3) = rc(3)
        end do

        ra(i) = -rc(3)+ac
        rv(i) = -rc(2)
        rd(i) = -rc(1)

        al = ac

    end do
    return

end subroutine rnmk

subroutine ranmk(acc,n,dt0,zeta,P,ra) bind(c)

    integer, intent(in) :: n
    real(8), intent(in) :: dt0, zeta
    real(8), intent(in) :: acc(n), P
    real(8), intent(out) :: ra(n)

    real(8), allocatable :: rv(:),rd(:)

    allocate(rv(n))
    allocate(rd(n))

    call rnmk(acc,n,dt0,zeta,P,ra,rv,rd)

    deallocate(rv)
    deallocate(rd)

    return
end subroutine ranmk

subroutine rfreq(acc,n,dt,zeta,P,ra,rv,rd) bind(c)

    integer, intent(in) :: n
    real(C_DOUBLE), intent(in) :: dt, zeta, P
    real(C_DOUBLE), intent(in) :: acc(n)
    real(C_DOUBLE), intent(out) :: ra(n),rv(n),rd(n)

    complex(C_DOUBLE_COMPLEX), allocatable :: af(:), raf(:), rvf(:), rdf(:)
    real(C_DOUBLE), allocatable :: a0(:), w(:)
    integer :: i, j, Nfft
    complex(C_DOUBLE_COMPLEX) :: IONE=cmplx(0.d0,1.d0)
    real(C_DOUBLE) :: w0i2, w0iwj, wj2, iNfft, w0
    type(C_PTR) :: plan, iplan

    Nfft = nextpow2(n)*4
    iNfft = 1.d0/dble(Nfft)

    allocate(a0(Nfft))
    allocate(af(Nfft))
    allocate(raf(Nfft))
    allocate(rvf(Nfft))
    allocate(rdf(Nfft))
    allocate(w(Nfft))

    call fftfreqs(Nfft,1.d0/dt,w)
    w = w*pi2
    w0 = pi2/P

    ra = 0.d0
    rv = 0.d0
    rd = 0.d0
    a0 = 0.d0
    a0(1:n) = acc

    plan  = fftw_plan_dft_r2c_1d(Nfft, a0, af, FFTW_ESTIMATE)
    iplan = fftw_plan_dft_c2r_1d(Nfft, af, ra, FFTW_ESTIMATE)
    call fftw_execute_dft_r2c(plan, a0, af)

    w0i2 = w0*w0
    do j = 1, Nfft, 1
        w0iwj = w0*w(j)
        wj2 = w(j)*w(j)
        raf(j) = af(j)*(w0i2+2.D0*zeta*w0iwj*IONE)/&
                (w0i2-wj2+2.d0*zeta*w0iwj*IONE)*iNfft
        rvf(j) = af(j)*(-w(j)*IONE/(w0i2-wj2+2.d0*zeta*w0iwj*IONE))*iNfft
        rdf(j) = af(j)*(-1.d0/(w0i2-wj2+2.d0*zeta*w0iwj*IONE))*iNfft
    end do
    call fftw_execute_dft_c2r(iplan, raf, ra)
    call fftw_execute_dft_c2r(iplan, rvf, rv)
    call fftw_execute_dft_c2r(iplan, rdf, rd)

    call fftw_destroy_plan(plan)
    call fftw_destroy_plan(iplan)

    deallocate( a0)
    deallocate( af)
    deallocate(raf)
    deallocate(rvf)
    deallocate(rdf)
    deallocate(  w)

    return

end subroutine rfreq

subroutine rafreq(acc,n,dt,zeta,P,ra) bind(c)

    integer, intent(in) :: n
    real(C_DOUBLE), intent(in) :: dt, zeta, P
    real(C_DOUBLE), intent(in) :: acc(n)
    real(C_DOUBLE), intent(out) :: ra(n)

    complex(C_DOUBLE_COMPLEX), allocatable :: af(:), raf(:)
    real(C_DOUBLE), allocatable :: a0(:), w(:)
    integer :: i, j, Nfft
    complex(C_DOUBLE_COMPLEX) :: IONE=cmplx(0.d0,1.d0)
    real(C_DOUBLE) :: w0i2, w0iwj, wj2, iNfft, w0
    type(C_PTR) :: plan, iplan

    Nfft = nextpow2(n)*4
    iNfft = 1.d0/dble(Nfft)

    allocate(a0(Nfft))
    allocate(af(Nfft))
    allocate(raf(Nfft))
    allocate(w(Nfft))

    call fftfreqs(Nfft,1.d0/dt,w)
    w = w*pi2
    w0 = pi2/P

    ra = 0.d0
    a0 = 0.d0
    a0(1:n) = acc

    plan  = fftw_plan_dft_r2c_1d(Nfft, a0, af, FFTW_ESTIMATE)
    iplan = fftw_plan_dft_c2r_1d(Nfft, af, ra, FFTW_ESTIMATE)
    call fftw_execute_dft_r2c(plan, a0, af)

    w0i2 = w0*w0
    do j = 1, Nfft, 1
        w0iwj = w0*w(j)
        wj2 = w(j)*w(j)
        raf(j) = af(j)*(w0i2+2.D0*zeta*w0iwj*IONE)/&
                (w0i2-wj2+2.d0*zeta*w0iwj*IONE)*iNfft
    end do
    call fftw_execute_dft_c2r(iplan, raf, ra)

    call fftw_destroy_plan(plan)
    call fftw_destroy_plan(iplan)

    deallocate( a0)
    deallocate( af)
    deallocate(raf)
    deallocate(  w)

    return

end subroutine rafreq

subroutine rmixed(acc,n,dt,zeta,P,ra,rv,rd) bind(c)

    integer, intent(in) :: n
    real(C_DOUBLE), intent(in) :: dt, zeta, P
    real(C_DOUBLE), intent(in) :: acc(n)
    real(C_DOUBLE), intent(out) :: ra(n),rv(n),rd(n)

    if (P<=MPR*dt) then
        call rfreq(acc,n,dt,zeta,P,ra,rv,rd)
    else
        call rnmk(acc,n,dt,zeta,P,ra,rv,rd)
    end if

end subroutine rmixed

subroutine ramixed(acc,n,dt,zeta,P,ra) bind(c)

    integer, intent(in) :: n
    real(C_DOUBLE), intent(in) :: dt, zeta, P
    real(C_DOUBLE), intent(in) :: acc(:)
    real(C_DOUBLE), intent(out) :: ra(:)

    if (P<=MPR*dt) then
        call rafreq(acc,n,dt,zeta,P,ra)
    else
        call ranmk(acc,n,dt,zeta,P,ra)
    end if

end subroutine ramixed

function rsimple(pe0,p0,phi,t,zeta) bind(c)
!     简谐激励下SDOF体系在t时刻的加速度响应
    real(C_DOUBLE), intent(in) :: pe0,p0,phi,t,zeta
    real(C_DOUBLE) :: rsimple

    real(C_DOUBLE) :: we, w
    real(C_DOUBLE) :: we2, w2, zeta2
    real(C_DOUBLE) :: we3, w3, zeta3
    real(C_DOUBLE) :: we4, w4, zeta4
    real(C_DOUBLE) :: sinphi, cosphi
    real(C_DOUBLE) :: two = 2.d0, three = 3.d0, four = 4.d0
    real(C_DOUBLE) :: eight = 8.d0, ten = 1.d1, twelve = 12.d0

    we = pi2/pe0
    w = pi2/p0
    we2 = we*we
    w2 = w*w
    we3 = we2*we
    w3 = w2*w
    we4 = we3*we
    w4 = w3*w

    zeta2 = zeta*zeta
    zeta3 = zeta2*zeta
    zeta4 = zeta3*zeta

    sinphi = sin(phi)
    cosphi = cos(phi)

    rsimple = cos(we*t+phi) - exp(-w*zeta*t)*((four*w*we3*zeta &
        -four*w*we3*zeta3)*exp(w*zeta*t)*sin(we*t+phi)+((two*we4 &
        -two*w2*we2)*zeta2-two*we4+two*w2*we2)*exp(w*zeta*t)*cos(we*t &
        +phi)+sqrt(four*w2-four*w2*zeta2)*(four*cosphi*w*we2*zeta3 &
        +two*sinphi*we3*zeta2+(cosphi*w3-three*cosphi*w*we2) &
        *zeta-sinphi*we3+sinphi*w2*we)*sin(sqrt(four*w2-four*w2 &
        *zeta2)*t/two)+(eight*cosphi*w2*we2*zeta4+four*sinphi*w &
        *we3*zeta3+(two*cosphi*w4-ten*cosphi*w2*we2)*zeta2-four &
        *sinphi*w*we3*zeta+two*cosphi*w2*we2-two*cosphi*w4) &
        *cos(sqrt(four*w2-four*w2*zeta2)*t/two))/(eight*w2*we2*zeta4 &
        +(two*we4-twelve*w2*we2+two*w4)*zeta2-two*we4+four*w2*we2 &
        -two*w4)

end function rsimple

subroutine rsimples(pe0,p0,phi,time,res,zeta,n) bind(c)
!     简谐激励下SDOF体系在time时刻的加速度响应
    integer(C_INT), intent(in) :: n
    real(C_DOUBLE), intent(in) :: pe0,p0,phi,zeta
    real(C_DOUBLE), intent(in) :: time(n)
    real(C_DOUBLE), intent(inout) :: res(n)
    real(C_DOUBLE) :: rsimple

    integer(C_INT) :: i

    real(C_DOUBLE) :: t
    real(C_DOUBLE) :: we, w
    real(C_DOUBLE) :: we2, w2, zeta2
    real(C_DOUBLE) :: we3, w3, zeta3
    real(C_DOUBLE) :: we4, w4, zeta4
    real(C_DOUBLE) :: sinphi, cosphi
    real(C_DOUBLE) :: two = 2.d0, three = 3.d0, four = 4.d0
    real(C_DOUBLE) :: eight = 8.d0, ten = 1.d1, twelve = 12.d0

    we = pi2/pe0
    w = pi2/p0
    we2 = we*we
    w2 = w*w
    we3 = we2*we
    w3 = w2*w
    we4 = we3*we
    w4 = w3*w

    zeta2 = zeta*zeta
    zeta3 = zeta2*zeta
    zeta4 = zeta3*zeta

    sinphi = sin(phi)
    cosphi = cos(phi)

    do i = 1, n, 1
        t = time(i)
        res(i) = cos(we*t+phi) - exp(-w*zeta*t)*((four*w*we3*zeta &
        -four*w*we3*zeta3)*exp(w*zeta*t)*sin(we*t+phi)+((two*we4 &
        -two*w2*we2)*zeta2-two*we4+two*w2*we2)*exp(w*zeta*t)*cos(we*t &
        +phi)+sqrt(four*w2-four*w2*zeta2)*(four*cosphi*w*we2*zeta3 &
        +two*sinphi*we3*zeta2+(cosphi*w3-three*cosphi*w*we2) &
        *zeta-sinphi*we3+sinphi*w2*we)*sin(sqrt(four*w2-four*w2 &
        *zeta2)*t/two)+(eight*cosphi*w2*we2*zeta4+four*sinphi*w &
        *we3*zeta3+(two*cosphi*w4-ten*cosphi*w2*we2)*zeta2-four &
        *sinphi*w*we3*zeta+two*cosphi*w2*we2-two*cosphi*w4) &
        *cos(sqrt(four*w2-four*w2*zeta2)*t/two))/(eight*w2*we2*zeta4 &
        +(two*we4-twelve*w2*we2+two*w4)*zeta2-two*we4+four*w2*we2 &
        -two*w4)
    end do

    return

end subroutine rsimples

subroutine rnmknl(acc,n,dt0,zeta,P,ra,rv,rd,para,np,model,tol,maxiter)  bind(c)
!
!  SUBROUTINE FOR COMPUTATION OF RESPONSE OF SDOF UNDER EARTHQUAKE
!  acc -----  GROUND ACC
!  N  -----  LENGTH OF acc
!  DEL  -----  ORIGINAL TIME INTERVAL
!  zeta  ----  DAMPING RATIO
!  PARA  ----  (/K0,FY,K1,XP,FP,KP,FEFF,KEFF/)
!     K0 -- initial stiffness
!     FY -- yielding force
!     K1 -- post-yielding stiffness
!     XP -- commited displacement of last step
!     FP -- commited force of last step
!     KP -- commited stiffness of last step
!   FEFF -- right-hand-side "force" of the equivalent static balance equation
!   KEFF -- left-hand-side "stiffness" of the equivalent static balance equation

    integer, intent(in) :: n
    real(8), intent(in) :: dt0, zeta, tol
    real(8), intent(in) :: acc(n), P
    real(8), intent(out) :: ra(n),rv(n),rd(n)
    integer, intent(in) :: np, maxiter, model
    real(8), intent(inout) :: para(np)

    real(8) :: rc(3),rl(3),ac,al,da,x
    real(8) :: k,c,w,feff,keff,beta,gamma,dt
    real(8) :: b1, b2, b3, b4, b5, b6, b7, b8
    integer :: i,j,r

    if ( dt0*MPR > P ) then
        r  = ceiling(MPR*dt0/P)
        dt = dt0/r
    else
        r = 1
        dt = dt0
    end if

    beta = 1.d0/4.d0
    gamma = 0.5d0

    b1 = 1.d0/(beta*dt*dt)
    b2 = 1.d0/(beta*dt)
    b3 = 1.d0/(beta*2.0d0)-1.d0
    b4 = gamma/(beta*dt)
    b5 = gamma/beta-1.d0
    b6 = 0.5d0*dt*(gamma/beta-2.0d0)
    b7 = dt*(1.0d0-gamma)
    b8 = dt*gamma

    w = 6.283185307179586d0/p
    k = w*w
    c = 2.0d0*zeta*w
    keff = b1+b4*c

    !
    !  INITIATION
    !
    para(1) = k
    para(8) = keff
    rl = [0.0d0,0.0d0,0.0d0]
    rc = [0.0d0,0.0d0,0.0d0]
    para(4:6) = (/rl(1),0.0d0,k/)
    al = 0.d0

    !
    !  COMPUTATION OF RESPONSE
    !
    do i = 1,n,1
        
        ac = acc(i)
        da = (ac - al)/dble(r)

        do j = 1, r, 1
            ac = al + da*j
    
            feff = ac+(b1*rl(1)+b2*rl(2)+b3*rl(3))&
            +c*(b4*rl(1)+b5*rl(2)+b6*rl(3))
            call solve_balance_equation(x,feff,para,np,model,tol,maxiter)

            rc(1) = x
            rc(2) = b4*(rc(1)-rl(1))-b5*rl(2)-b6*rl(3)
            rc(3) = ac-para(5)-c*rc(2)

            rl(1) = rc(1)
            rl(2) = rc(2)
            rl(3) = rc(3)
        end do

        ra(i) = -rc(3)+ac
        rv(i) = -rc(2)
        rd(i) = -rc(1)

        al = ac
    end do
    return

end subroutine rnmknl

subroutine newmark_nl(acc,n,dt0,zeta,P,SPA,SPI,SPV,SPD,SPE,para,np,model,tol,maxiter) bind(c)
!
!  SUBROUTINE FOR COMPUTATION OF RESPONSE OF SDOF UNDER EARTHQUAKE
!  acc -----  GROUND ACC
!  N  -----  LENGTH OF acc
!  DEL  -----  ORIGINAL TIME INTERVAL
!  zeta  ----  DAMPING RATIO
!  PARA  ----  (/K0,FY,K1,XP,FP,KP,FEFF,KEFF/)
!     K0 -- initial stiffness
!     FY -- yielding force
!     K1 -- post-yielding stiffness
!     XP -- commited displacement of last step
!     FP -- commited force of last step
!     KP -- commited stiffness of last step
!   FEFF -- right-hand-side "force" of the equivalent static balance equation
!   KEFF -- left-hand-side "stiffness" of the equivalent static balance equation

    integer, intent(in) :: n
    real(8), intent(in) :: dt0, zeta, tol
    real(8), intent(in) :: acc(n), P
    real(8), intent(out) :: SPA,SPV,SPD,SPE
    integer, intent(out) :: SPI
    integer, intent(in) :: np, maxiter, model
    real(8), intent(inout) :: para(np)

    real(8) :: rc(3),rl(3),ac,al,da,x,aa,rvl,ein
    real(8) :: k,c,w,feff,keff,beta,gamma,dt
    real(8) :: b1, b2, b3, b4, b5, b6, b7, b8
    integer :: i,j,r


    if ( dt0*MPR > P ) then
        r  = ceiling(MPR*dt0/P)
        dt = dt0/r
    else
        r = 1
        dt = dt0
    end if

    beta = 1.d0/4.d0
    gamma = 0.5d0

    b1 = 1.d0/(beta*dt*dt)
    b2 = 1.d0/(beta*dt)
    b3 = 1.d0/(beta*2.0d0)-1.d0
    b4 = gamma/(beta*dt)
    b5 = gamma/beta-1.d0
    b6 = 0.5d0*dt*(gamma/beta-2.0d0)
    b7 = dt*(1.0d0-gamma)
    b8 = dt*gamma

    w = 6.283185307179586d0/p
    k = w*w
    c = 2.0d0*zeta*w
    keff = b1+b4*c

    !
    !  INITIATION
    !
    para(1) = k
    para(8) = keff
    rl = [0.0d0,0.0d0,0.0d0]
    rc = [0.0d0,0.0d0,0.0d0]
    para(4:6) = (/rl(1),0.0d0,k/)
    al = 0.d0
    rvl = 0.d0

    SPA = 0.d0
    SPV = 0.d0
    SPD = 0.d0
    SPE = 0.d0
    SPI = 0

    !
    !  COMPUTATION OF RESPONSE
    !
    do i = 1,n,1
        
        ac = acc(i)
        da = (ac - al)/dble(r)

        do j = 1, r, 1
            ac = al + da*j
    
            feff = ac+(b1*rl(1)+b2*rl(2)+b3*rl(3))&
            +c*(b4*rl(1)+b5*rl(2)+b6*rl(3))
            call solve_balance_equation(x,feff,para,np,model,tol,maxiter)

            rc(1) = x
            rc(2) = b4*(rc(1)-rl(1))-b5*rl(2)-b6*rl(3)
            rc(3) = ac-para(5)-c*rc(2)

            rl(1) = rc(1)
            rl(2) = rc(2)
            rl(3) = rc(3)
        end do

        SPE = SPE+0.5d0*dt0*(rvl*al+rc(2)*ac)
        aa = -rc(3)+ac

        if (abs(aa) > abs(SPA)) then
            SPA=aa
            SPI=i
        endif

        if (abs(rc(2)) > abs(SPV)) then
            SPV=rc(2)
        endif

        if (abs(rc(1)) > abs(SPD)) then
            SPD=rc(1)
        endif

        al = ac
        rvl = rc(2)

    end do

    SPE = sqrt(2.d0*SPE)

    return

end subroutine newmark_nl

subroutine solve_balance_equation(x,feff,para,np,model,tol,maxiter) bind(c)

    integer, intent(in) :: np, maxiter, model
    real(8), intent(in) :: feff, tol
    real(8), intent(inout) :: para(np)
    real(8), intent(out) :: x

    real(8) :: x0

    !     估计变形值
    x0 = para(4)+(feff-para(7))/(para(6)+para(8))
    !     施加荷载
    para(7) = feff
    !     迭代求解变形
    call newton(x,x0,tol,maxiter,para,np,model)
    !     更新状态变量(para(4:6))
    call update_para(x,para,np,model)

    return

end subroutine solve_balance_equation

subroutine newton(x,x0,tol,maxiter,para,np,model) bind(c)
!     牛顿-拉普森法求解非线性方程
!     输入参数：
!           x0 : 解的初始估计值
!           tol : 收敛容差(相对)
!           maxiter : 最大迭代次数
!           para : 一维数组, 包含非线性参数及历史状态变量, len(para)=np
!     输出：
!           x : 方程的解

    integer, intent(in) :: maxiter,np,model
    real(8), intent(out) :: x0,x
    real(8), intent(in) :: tol
    real(8), intent(inout) :: para(np)
    
    real(8) :: fval, kc
    integer :: iter
    !
    !     external func,f_prime
    !     newton-rapheson method
    !f2py intent(in) x0,tol,maxiter,para,np
    !f2py intent(out) x
    kc = 0.0d0
    fval = 0.0d0
    x = x0
    do iter = 1,maxiter,1

        call balance_func(fval,kc,x,para,np,model)

        if (kc .ne. 0.d0) then
            x = x0 - fval/kc
            if (abs(x - x0) .lt. tol) then
                exit
            endif
            x0 = x
        else
            write(*,*) "Derivative was zero. Stop the Iteration."
            x = x0
            exit
        endif

    end do
    !       write(*,*) "iter = ",iter,"kc = ",kc

end subroutine newton

subroutine update_para(x,para,np,model) bind(c)
!     更新历史状态变量

    integer, intent(in) :: np, model
    real(8), intent(in) :: x
    real(8), intent(inout) :: para(np)

    real(8) :: fval, kc

    fval = 0.0d0
    kc = 0.0d0
    call balance_func(fval,kc,x,para,np,model)

    para(4) = x
    para(5) = fval-para(8)*x+para(7)
    para(6) = kc-para(8)

end subroutine update_para

subroutine balance_func(fval,kc,x,para,np,model) bind(c)
!     已知X, 计算方程F(X)=0中F(X)的值及F'(X)的值KC
!     输入参数：
!           X ---- 变形值
!           PARA,NP ---- 一维数组, 包含非线性参数及历史状态变量,
!                       LEN(PARA)=NP
!           MODEL ---- 非线性滞回模型
!                      0 ---- 双线性(DEFAULT)
!                      1 ---- CLOUGH 退化双线性
!                      2 ---- 武田 退化三线性
!                      3 ---- 无残余变形双线性(卸载指向原点)
!     输入结果：
!           FVAL ---- 平衡力残差
!           KC ---- 当前变形对应的刚度值(切线斜率)
!
    integer, intent(in) :: np, model
    real(8), intent(in) :: x
    real(8), intent(in) :: para(np)

    real(8), intent(out) :: fval, kc

    real(8) :: k0,fy,k1,xp,fp,kp,feff,ck,dx,f,f_try,bup,bdn

    k0 = para(1)
    fy = para(2)
    k1 = para(3)

    xp = para(4)
    fp = para(5)
    kp = para(6)

    feff = para(7)
    ck = para(8)

    dx = x-xp
    f_try = fp+dx*kp

    select case (model)
        case (1)
        case (2)
        case (3)
        case default
            bup = fy+k1*(x-fy/k0)
            bdn = -fy+k1*(x+fy/k0)

            if (dx .gt. 0.0d0) then
                if (f_try .lt. bup) then
                    kc = k0
                    f = min(fp+dx*kc,bup)
                else
                    f = bup
                    kc = k1
                end if

            elseif (dx .lt. 0.0d0) then
                if (f_try .gt. bdn) then
                    kc = k0
                    f = max(fp+dx*kc,bdn)
                else
                    f = bdn
                    kc = k1
                end if
            else
                kc = kp
                f = fp
            end if
    end select

    kc = kc+ck
    fval = f+ck*x-feff

end subroutine balance_func

subroutine energy(acc,n,dt,zeta,P,ra,rv,rd,Ek,Ez,Es,Eh,Ein,BL) bind(c)
    
    integer, intent(in) :: n
    real(8), intent(in) :: acc(n), ra(n), rv(n), rd(n)
    real(8), intent(in) :: dt, zeta, P
    real(8), intent(out) :: Ek(n),Ez(n),Es(n),Eh(n),Ein(n),BL(n)

    integer :: i
    real(8) :: w, fs(2), Ea, c

    w = pi2/P
    c = 2.d0*zeta*w

    Ek(1) = rv(1)*rv(1)*0.5D0
    Ez(1) = 0.5D0*dt*(0.d0+c*rv(1)*rv(1))
    fs = [0.D0, -ra(1)-c*rv(1)]
    Ea = 0.5D0*dt*(0.d0+fs(2)*rv(1))
    Es(1) = fs(2)*fs(2)*0.5D0/(w*w)
    Eh(1) = Ea - Es(1)
    Ein(1) = 0.5D0*dt*(0.d0-acc(1)*rv(1))
    BL(1) = Ein(1) - Ek(1) - Ez(1) - Ea

    do i = 2, n, 1
        Ek(i) = rv(i)*rv(i)*0.5D0
        Ez(i) = Ez(i-1) + 0.5D0*dt*c*(rv(i-1)*rv(i-1)+rv(i)*rv(i))

        fs(1) = - ra(i-1) - c*rv(i-1)
        fs(2) = - ra(i) - c*rv(i)
        Ea = Ea + 0.5D0*dt*(fs(1)*rv(i-1)+fs(2)*rv(i))
        Es(i) = fs(2)*fs(2)*0.5D0/(w*w)
        Eh(i) = Ea - Es(i)

        Ein(i) = Ein(i-1) - 0.5D0*dt*(acc(i-1)*rv(i-1)+acc(i)*rv(i))
        BL(i) = Ein(i) - Ek(i) - Ez(i) - Ea

    end do
    
end subroutine energy

end module response
